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Risks and damages associated with lava flows propagation (for 
instance the most recent Etna eruptions) require a quantitative de- 
scription of this phenomenon and a reliable forecasting of lava flow 
paths. Due to the high complexity of these processes, numerical 
solution of the complete conservation equations for real lava flows 
is often practically impossible. To overcome the computational 
difficulties, simplified models are usually adopted, including 1-D 
models and cellular automata. In this work we propose a simplified 
2D model based on the conservation equations for lava thickness 
and depth-averaged velocities and temperature which result in first 
order partial differential equations. The proposed approach repre- 
sents a good compromise between the full 3-D description and the 
need to decrease the computational time. The method was satisfac- 
torily applied to reproduce some analytical solutions and to simulate 
a real lava flow event occurred during the 1991-93 Etna eruption. 
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where h is the fluid depth measured from the altitude of the terrain 
surface H (bed), (17, V) = 1/h J^ +h u ( x, y, z)dz are the depth- 
averaged fluid velocity components, j3ij are correction factors (in 
the range 0.5-1.5) and 7 is a dimensionless friction coefficient de- 
pending on the fluid rheology and on the properties of both flow and 
bed. The gradients dH/dxi indicate the channel bottom slopes in 
both directions x and y (Xi — x,y). The terms on the right sides 
represent the so-called source terms. 

In the case of lava, the viscosity is strongly temperature dependent. 
For this reason, besides the equations (1), (2) and (3), it is necessary 
to solve the equation for the energy conservation. From a com- 
putational point of view, the temperature equation is similar to the 
pollutant transport equation [Monthe et al, 1999; LeVeque, 2002], 
We propose the following heuristic equation for the depth-averaged 
temperature T(x, y) — 1/h J^ +h T(x, y, z)dz: 



1. Introduction 

Depth averaged flow models based on the so-called shallow wa- 
ter equations (SWE) were firstly introduced by De Saint Venant in 
1864 and Boussinesq in 1872. Nowdays, applications of the shal- 
low water equations include a wide range of problems which have 
important implications for hazard assessment, from flood simula- 
tion [Burguete et al. , 2002] to tsunamis propagation [Heinrich et al, , 
2001]. 

In this paper we propose a generalized set of depth averaged equa- 
tions, including an energy equation, to describe lava flow propaga- 
tion. We considered lava flow as channelized, i.e. moving lava has 
a non-continuous roof and the top represents a free surface open to 
the atmosphere. 
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where T c and T env are the temperatures of the lava-ground inter- 
face and of the external environment respectively, and f3n, £, W, 
Ti and K. are a set of semi-empirical parameters. Terms on the 
right side of the equation (4) represent the radiative, convective and 
conductive exchanges respectively, while the last term is due to the 
viscous heating. Moreover, a simple exponential relationship be- 
tween magma viscosity and temperature was assumed [Costa and 
Macedonio, 2002]: 



2. Model description 

The model is based on depth-averaged equations obtained by 
integrating mass, momentum and energy equations over the fluid 
depth, from the bottom up to the free surface. This approach is 
valid in the limit H 2 /L 2 <C 1 (where H* is the undisturbed fluid 
height and L» the characteristic wave length scale in the flow direc- 
tion). This means that we are dealing with very long waves or with 
"shallow water". 

Assuming an incompressible homogeneous fluid and a hydrostatic 
pressure distribution, the shallow water equations for an uniform or 
gradually varied flow are given by: 
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where b is an appropriate rheological parameter and fj, r is the viscos- 
ity value at the reference temperature T r (for instance, T r = To with 
To equal to the emission temperature at the vent). For the descrip- 
tion of a thermal balance in lava flows, similar to the equation (4) 
see Keszthely and Self [1998]. We do not explicitly accounted 
for crystallization and crystallinity-dependence of the viscosity, but 
they are implicitly considered in the determination of the rheologi- 
cal parameters in (5). Concerning the coefficient 7 which appears 
in the equations (2) and (3), we propose a relationship similar to 
that used in the viscous regime [Gerbeau and Perthame, 2001; Fer- 
rari and Saleri, 2004]: 7 = k*/[1 + K t h/(3v r )], where k* is the 
Navier friction coefficient, v r — fj, r /p and p — fluid density. This 
relationship permits in principle to consider different and general 
wall friction conditions and, for instance, the possibility to include 
viscous heating effects on lava flow velocity [Costa and Macedo- 
nio, 2003] by choosing the appropriate «» parameterization. By 
considering the viscosity dependence on temperature(5) and, for 
simplicity, the limit n,h/ (3v r ) S> 1, we obtain: 
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In the following, we estimate the other parameters introduced in (4) 
evaluating the corresponding terms of the complete averaged energy 
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equation. The heat transfer coefficient H is roughly estimated from 
the term k J^ +h V 2 T(x,y, z)dz: 



H w nn/h 



(7) 



where tt = k/(pc p ) is the thermal diffusivity (k is the thermal con- 
ductivity and c p the specific heat) and we approximated the char- 
acteristic thermal boundary layer length as a fraction of the total 
thickness: St ~ h/n where n depends on the temperature profile 
(n ~ 4 ~ ^JvJk). 

According to Fieri andBaloga [1986]'s study, for the radiative term, 
we assumed: 
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where e is the emissivity, a the Stephan-Boltzmann constant (a = 
5.67 • 10~ 8 WmT 2 K~ 4 ) and / is the fractional area of the exposed 
inner core [Crisp and Baloga, 1990]. For simplicity, in this version 
of the model we assumed / as a constant. In real lava flows / may 
change with time and space / = /(x,t) and, in principle, it can 
be estimated from field measurements or remote sensing. Further 
studies should investigate the sensivity of the model with the tem- 
poral and spatial changes of this quantity. 
For the convective term, we adopted [Keszthely and Self, 1998]: 

W « Xf/(pc p ) (9) 

where A is the atmospheric heat transfer coefficient. 

Finally, for the viscous heating term, we approximate the order of 

magnitude of the quantity $ = l/(pc p ) j^ +h fi(dv/dz) 2 dz as 

p, r e~ b( - T ~ Tr \U 2 + V 2 )m/h, where we approximated the charac- 
teristic velocity boundary layer as S v « h/m; hence: 

K, « mpr/{pCph) (10) 

where in the case of a parabolic velocity profile m = 12 [Shah and 
Pearson, 1974]. 

By using the approximations and parameterizations described 
above, we obtain the final system of equations we solve by means 
of the numerical method described in the Section 3. 



Figure 1. Longitudinal profiles of the channel center velocity 
and temperature, at t = 2500 s. Dashed and continuous lines 
indicate analytical and numerical results, respectively. Channel 
dimensions: 50 m wide, 1000 m long; Slope: 0.1, T env = K; 
T = 1353 K, Flow rate: Q = 12.5 m 3 /s; Ax = Ay = 5 m. 
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3. The numerical method 

The numerical solution of the equations (1), (2), (3) 
and (4), was achieved by using an algorithm based on 
the software package CLAWPACK (available on the web 
at http : //www . amath . Washington . edu/ ~ rl j / claw- 
pack . html). CLAWPACK is a public domain software package 
designed to compute numerical solutions to hyperbolic partial dif- 
ferential equations using a wave propagation approach [LeVeque, 
2002]. 

The CLAWPACK routines were generalized in order to treat the 
viscous friction source term and to solve the energy equation (4). 
The modelling of lava flow over an initially dry downstream region 
(dry bed problem) was approached following the method described 
in Monthe et al. [1999]. All the source terms in the governing equa- 
tions were treated using a Godunov splitting method and, since as 
a simple explicit discretization leads to numerical instabilities [e.g. 
Ambrosi, 1999; Monthe et al, 1999], all terms were discretized us- 
ing a semi-implicit scheme. For instance, the source term in the 
equation (2) was discretized as below: 

n n + l ~ In _ , 6H 3lS r qn + l _(,(T n -T r ) 

At ~ ~ 9 n ~dx~ hl~ e 

where pedice n indicates the quantities at the time t„, and q n — 
U n h n . The other source terms were discretized by using a similar 
approach. 

Before the application, the algorithm was tested by simulating 
some cases for which analytical solutions are known. In fact, con- 
sidering the flow of a quasi-unconfined layer of viscous liquid on 



Figure 2. Longitudinal temperature profiles at t — 2500 s ob- 
tained using the same parameters of Figure 1 and T env = 300 K, 
T c = 1173 K, / = 0.5, e = 0.8, n = 4, m = 12. VH. = Vis- 
cous Heating. 



an inclined plane, with the energy and the momentum equations 
decoupled (i.e. with b — K~ ) and in the steady state limit, the 
equations (1), (2), (3) and (4) admit the following analytical rela- 
tionships [Keszthely and Self , 1998; Fieri and Baloga, 1986]: 

1i = -g?3sina/(3^ r ) q 3 = qilT^ 3 + 3£(y - yo)/q 2 r W ' i 

(11) 

where qi — h, q? = hV, qs = hT, a is the channel slope and 
(y — yo) represents the distance from the vent. Figure 1 shows 
the comparison between the analytical and numerical relationships. 
Simulation results have shown a good agreement with an error less 
than 1% for the conservative variables h, hV and hT and, within 
a few % for the non-conservative variable V and T. Moreover, in 
order to estimate the importance of each term on the right side of the 
equation (4), we considered the same geometry of the simple slope 
flow as above and the typical values reported in the caption of the 
Figure 2. Results, plotted in the Figure 2, show that radiative cooling 
is the main heat loss mechanism, while conductive and atmospheric 
convective cooling is less important but, for the parameter values 
used here, conductive loss is comparable with convection cooling. 
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its typical values since, for numerical reasons, we need to limit the 
maximum viscosity value. 

The parameters reported in Table 1 give the following typical values: 
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Figure 3. Longitudinal thickness profiles at t — 1200 s obtained 
considering the same parameters reported in Figure 2. 



Viscous heating effect can be neglected in terms of mean lava tem- 
perature (in the simulated case it produces a increase of a few °C 
for a distance of 1 km), although, in certain conditions, it could 
be more important and determinant in the choosing the appropriate 
wall conditions and exchange coefficients for both momentum and 
energy [Costa andMacedonio, 2003]. About effects of the coupling 
between momentum and energy equations, we can see a non-zero b 
is important to determine the longitudinal variation of the lava flow 
thickness (see Figure 3), although it increases slightly the cooling 
beyond certain distances. Figure 3 shows as the velocity decrease 
due to the longitudinal viscosity increase is able to cause a longitu- 
dinal rise of the lava thickness because of the viscosity temperature 
dependence. 

4. Application to Etna lava flows 

In this section, as an application, we reported simulation results 
of the initial phases of the 1991-1993 Etna eruption for which some 
field data for input and comparison are available [Calvari et al, 
1994]. In particular we simulated the second phase occurred from 
the 3 rd up to the 10 th January 1992. In order to estimate previously 
introduced semi-empirical parameters, we considered the typical 
magma parameters reported in Table 1 partially derived from data 
of Calvari et al. [1994]. We assumed as representative an effec- 
tive viscosity of 10 3 Pa s at an estimated vent temperature of about 
1353 K and b ta 0.02 KT 1 that, for a cooling of about 100 K, re- 
produces the observed viscosities of the order of 10 4 Pas [Calvari 
et al, 1994]. Other parameters were chosen within typical ranges: 
/ = 0.1 (between 0.01 and 1 [Keszthely and Self, 1998]), and 
e = 0.8 (between 0.6 and 0.9 [Neri, 1998]). T c is set higher than 




H~3/hx 10 _6 ms _1 
£ « 1.5 x 10 _15 ms _:l K -3 
Wsi2x 10 _6 ms _1 



(12) 



K ~ 4/7i x 10" 



where, for our aim in this application, we set T env = 300 K, n = 4, 
m = 12 and /3y = 1. 



Table 1. Parameters characteristic of Etna lava. 
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Figure 4. Simulated lava thickness of the 3 rd and 4 th January 
Etna lava flow. 



As topographic basis, we used the digital data files of the Etna 
maps with a 1 : 10000 scale available at the Osservatorio Vesuviano- 
INGV web site at http : / /venus . ov . ingv . it (the used spa- 
tial grid resolution was Aa; = Ay = 25 m). For the second phase, 
we considered an ephemeral vent sited in Piano del Trifoglietto at 
the UTM coordinates (503795; 4174843). Finally, for the period 
3-10 January 1992, we considered a constant average lava flow rate 
of 16 m 3 /s (ranging from 8 to 25 m 3 /s) [Calvari etal., 1994; Barberi 
etal, 1993]. 

The first phase of the eruption corresponded with the initial spread- 
ing of the lava flows on Piano del Trifoglietto. On the 3 rd January 
1992 a new lava flow that overlapped the older lava lows, became 
an independent branch. By the evening it covered more than 1 km. 
The day after the front reached Mt. Calanna. One branch continued 
to move to the south of Mt. Calanna and one branch turned to the 
north then to the east (see Figure 3 of Calvari et al. [1994]). Be- 
cause of a significantly decrease of lava supply, the southern lava 
flow stopped in Val Calanna. On January 7* the northern lava lobe 
touched the southern one and then merged [Calvari et al, 1994]. 
In Figure 4 the simulated lava flow at the end of the second phase 
is shown. The model is able to reproduce semi-quantitatively the 
behaviour of the real lava flow and the order of magnitude of the 
quantities involved such as thickness, temperature and the time of 
front propagation of the lava flow. Although we introduced different 
simplifications and we considered an arduous case encompassing 
both a large viscous friction term and complex rough topography, 
the simulation and real lava flows show strikingly similar dynamics 
and thermal pattern evolution. Nevertheless the model presented 
in this paper remains an initial model of lava flow emplacement 
using SWE. Future improvements are expected by refining the the 
computational performance of the model and the formulation of the 
parameters. 

5. Limitations 

This methodology is based on vertical averages and therefore it 
cannot be rigorously valid for every conceivable application. We 
stress that the model is based on the basic assumptions of (1) small 
vertical scale relative to horizontal (//»/L» <C 1), (2) homogeneous 
incompressible fluid, (3) hydrostatic pressure distribution, (4) slow 
vertical variations. 
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Concerning the computational method, the principal limit is related 
to the numerical treatment we used here for the source terms arising 
from topography and viscous friction. In particular since the ac- 
tual topographies may contain abrupt variations, the slope term that 
appears in the equations (2) and (3) can become infinite in correspon- 
dence of discontinuities leading to numerical oscillations, diffusion, 
smearing and non-physical solutions [LeVeque, 1998; Alcrudo and 
Benkhaldoun, 2001; Chinnayya etal, 2004]. Also the friction term 
must be carefully treated. In fact, if the characteristic time of the 
source term is much smaller than the characteristic time of the con- 
vective part of the equations, the problem is said to be stiff and the 
classical splitting method may provide erroneous physical solutions 
on coarse meshes [Chinnayya et al. , 2004] . To avoid these problems 
a trivial solution is using a very small time step, which results in 
long computational times. In the next version of the model, this 
limit could be overcome by applying directly a method based on the 
solution of the inhomogeneous Riemann problem with source term 
instead of applying the splitting method [Chinnayya et al, 2004; 
George, 2004]. 

6. Conclusion 

A new general computational model for lava flow propagation 
based on the solution of depth-averaged equations for mass, momen- 
tum and energy equation was described. This approach appears to 
be a robust physical description and a good compromise between 
the full 3-D simulation and the necessity to decrease the compu- 
tational time. The model was satisfactorily applied to reproduce 
some analytical solutions and to simulate a real lava flow event oc- 
curred during the 1991-93 Etna eruption. The good performance 
obtained in this preliminary version of the model makes this ap- 
proach a potential tool to forecast reliably lava flow paths to use for 
risk mitigation, although the used algorithm should be improved for 
a better treatment of the source terms. 
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